#include "plant_SPDE2D.h"

// --------------------------------------------------------------
plant_SPDE2D::plant_SPDE2D(int M, int L,double idt){
	int i,m,n;
	mode = M;	
	int idim[2] = {1,1}; // Input#1 = T(K), Input#2 = W (layer/s)
	int odim[4];
	odim[0] = (mode+1)*(mode+1); // The mode with sin*sin
	odim[1] = (mode+1)*(mode+1); // 
	odim[2] = (mode+1)*(mode+1);
	odim[3] = (mode+1)*(mode+1);
	block::initialize(2,idim,4,odim,idt);

	LatticeSize = L;

//	z1 = block::outputs_next[0];
//	z2 = block::outputs_next[1];
//	z3 = block::outputs_next[2];
//	z4 = block::outputs_next[3];

	z1 = new double[(mode+1)*(mode+1)];
	z2 = new double[(mode+1)*(mode+1)];
	z3 = new double[(mode+1)*(mode+1)];
	z4 = new double[(mode+1)*(mode+1)];

	lambda = new double*[mode+1];
	for(m=0;m<mode+1;m++){
		lambda[m] = new double[mode+1];
	}

	for(i=0;i<4;i++){
		Km2[i] = new double[(mode+1)*(mode+1)];
	}

//	double dphi;
//	int ix,iy,ix_next;
	int index;
	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			index = index2D(m,n);
			if(m==0||2*m==LatticeSize||n==0||2*n==LatticeSize){
				Km2[0][index] = 0.0;
			}else{
				Km2[0][index] = 4.0/(pi*pi)*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize);
			}

			if(m==0){
				Km2[1][index] = 0.0;
			}else{
				Km2[1][index] = 4.0/(pi*pi)*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize);
			}

			if(m==0||2*m==LatticeSize){
				Km2[2][index] = 0.0;
			}else{
				Km2[2][index] = 4.0/(pi*pi)*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize);
			}

			if(n==0 || 2*n==LatticeSize || m==0){
				Km2[3][index] = 0.0;
			}else{
				Km2[3][index] = 4.0/(pi*pi)*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize);
			}
		}
	}
/*	ofstream fout_Km2;
	fout_Km2.open("Output_Km2.dat");
	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			fout_Km2 << Km2[1][index2D(m,n)] << " ";
		}
		fout_Km2 << "\n";
	}
	fout_Km2.close();*/

	dt_conservative_factor = 0.1;
	sqrt_dt = sqrt(block::dt);
	
	// Initialize random noize generator
	gsl_rng_env_setup();
	noise_type = gsl_rng_default;
	noise = gsl_rng_alloc(noise_type);
	gsl_rng_set(noise,time_seed());
	
	ModelC = new model_C_001();
	ModelC2 = new model_C2_004();
	ModelSigma2 = new model_sigma2_004();

	reset();
}
// --------------------------------------------------------------
plant_SPDE2D::~plant_SPDE2D(){
	int m;
	delete []z1;
	delete []z2;
	delete []z3;
	delete []z4;
	for(int i=0;i<4;i++)
	{
		delete []Km2[i];
	}

	for(m=0;m<mode+1;m++){
		delete []lambda[m];
	}
	delete []lambda;
	gsl_rng_free(noise);

	delete ModelC;
	delete ModelC2;
	delete ModelSigma2;
}
// --------------------------------------------------------------
void plant_SPDE2D::reset(){
	block::time = dt;
	int m,n;
	int index;
	for(m=0;m<mode+1;m++){
		for(n=0;n<mode+1;n++){
			index = index2D(m,n);
			z1[index] = 0.0;
			z2[index] = 0.0;
			z3[index] = 0.0;
			z4[index] = 0.0;
			outputs[0][index] = 0.0;
			outputs[1][index] = 0.0;
			outputs[2][index] = 0.0;
			outputs[3][index] = 0.0;
			outputs_next[0][index] = 0.0;
			outputs_next[1][index] = 0.0;
			outputs_next[2][index] = 0.0;
			outputs_next[3][index] = 0.0;
		}
	}
	
	T = 610;
	W = 0.5;
	double T_W[2];
	// c = eq_c(T,W);
	// c2 =eq_c2(T,W);
	// sigma = sqrt(eq_sigma2(T,W));
	T_W[0] = T;
	T_W[1] = W;
	c = ModelC->value(T_W);
	c2= ModelC2->value(T_W);
	sigma = sqrt(ModelSigma2->value(T_W));

//	c = 0.5;
//	c2 = 0.1;
//	sigma = 0.1;

	for(m=0;m<mode+1;m++){
		for(n=0;n<mode+1;n++){
			lambda[m][n] = -4*c2*(m*m+n*n);
		}
	}
}
// --------------------------------------------------------------
void plant_SPDE2D::update(double sysTime){
	double T = inputs[0][0];
	double W = inputs[1][0];
	update(sysTime,T,W);
}
// --------------------------------------------------------------
void plant_SPDE2D::update(double sysTime, double iT, double iW){	
	if(time>sysTime+dt) return;

	int m,n;
	int index;
	double c_200;
//	double dt_con;
//	double sqrt_dt_con;
//	double time_mn;
	double temp1;
	double temp2;
	
	if(time<=sysTime+dt){
		block::output();
	}
	while(time<=sysTime+dt){		
		if(iT!=T || iW!=W){
			T = iT;
			W = iW;
			//c = eq_c(T,W);
			//c2 = eq_c2(T,W);
			double T_W[2];
			T_W[0] = T;
			T_W[1] = W;
			c = ModelC->value(T_W);
			c2 = ModelC2->value(T_W);			
			//sigma = sqrt(eq_sigma2(T,W));
			sigma = sqrt(ModelSigma2->value(T_W));
			for(m=0;m<mode+1;m++){
				for(n=0;n<mode+1;n++){
					lambda[m][n] = -4*c2*(m*m+n*n);
				}
			}
		}
		
		c_200 = LatticeSize*c;
		for(m=0;m<=mode;m++){
			for(n=0;n<=mode;n++){
				index = index2D(m,n);
				temp1 = exp(lambda[m][n]*dt);
				temp2 = sigma*sqrt((exp(2.0*lambda[m][n]*dt)-1.0)/(2.0*lambda[m][n]));
				if(m>0&&n>0){
					z1[index] = z1[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
				}

				if(m==0&&n==0){
					z2[index] = z2[index]+c_200*dt+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma*sqrt_dt;
				}
				else{
					z2[index] = z2[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
				}

				if(m>0){
					z3[index] = z3[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
				}

				if(n>0){
					z4[index] = z4[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
				}
			}
		}			
		time += dt;
	}

	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			index = index2D(m,n);
			block::outputs_next[0][index] = z1[index]*z1[index];
			block::outputs_next[1][index] = z2[index]*z2[index];
			block::outputs_next[2][index] = z3[index]*z3[index];
			block::outputs_next[3][index] = z4[index]*z4[index];
		}
	}
	return;
}
// --------------------------------------------------------------
void plant_SPDE2D::output(string FileName){
	ofstream fout;	
	fout.open(FileName.c_str());
	fout << "Mode = " << mode;
	int i;

	int dim = (mode+1)*(mode+1);

	for(i=0;i<dim;i++){
		fout << z1[i] << " ";
	}
	fout << "\n";

	for(i=0;i<dim;i++){
		fout << z2[i] << " ";
	}
	fout << "\n";

	for(i=0;i<dim;i++){
		fout << z3[i] << " ";
	}
	fout << "\n";

	for(i=0;i<dim;i++){
		fout << z4[i] << " ";
	}
}
// --------------------------------------------------------------
double plant_SPDE2D::getM2(){
	int m,n,index;
	double m2 = 0.0;
	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			index = index2D(m,n);
			m2 +=  Km2[0][index]*outputs[0][index]
			      +Km2[1][index]*outputs[1][index]
			      +Km2[2][index]*outputs[2][index]
            	  +Km2[3][index]*outputs[3][index];
		}
	}
	return m2;
}
// --------------------------------------------------------------
double plant_SPDE2D::getR2(){
	int m,n,index;
	double r2 = 0.0;
	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			if(m!=0||n!=0){
				index = index2D(m,n);
				r2 +=  outputs[0][index]
				      +outputs[1][index]
				      +outputs[2][index]
				      +outputs[3][index];
				//r2 += outputs[1][index];
			}
		}
	}
	//r2 = outputs[1][index2D(0,1)];
	r2 /= (pi*pi);
	return r2;
}
// --------------------------------------------------------------
/*double plant_SPDE2D::eq_c(double T,double W)
{
	return W*(1-k_w*pow(W,-a_w)*exp(k_B*T/E_w));
}

double plant_SPDE2D::eq_c2(double T,double W)
{
	return k_c*exp(k_B*T/E_c)*pow(W,-a_c)/(LatticeSize*LatticeSize);
}

double plant_SPDE2D::eq_sigma2(double T,double W)
{
//	double index = (a_t+k_t*W)*T-a_v-k_v*W;
//	return pi*pi/(LatticeSize*LatticeSize)*(1+exp(index));

	return (T0_W0+T0_W1*W)+(T1_W0+T1_W1*W)*T+(T2_W0+T2_W1*W)*T*T+(T3_W0+T3_W1*W)*T*T*T;
}*/
// --------------------------------------------------------------
unsigned plant_SPDE2D::time_seed (void)
{
  time_t timeval;       /* Current time. */
  unsigned char *ptr;   /* Type punned pointed into timeval. */
  unsigned seed;        /* Generated seed. */
  size_t i;

  timeval = std::time (NULL);
  ptr = (unsigned char *) &timeval;

  seed = 0;
  for (i = 0; i < sizeof timeval; i++)
    seed = seed * (0xff + 2u) + ptr[i];

  return seed;
}
